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ABSTRACT 

Information theory can develop a technique which takes experimentally deter- 
mined numbers and produces a uniquely specified "lies*' density model satisfying 
those numbers. This technique does not depend on previous density models for 
starting conditions; it is self-starting. A model was generated using five nu- 
merical parameters: the mass of the earth, its moment of inertia, three zero- 
node torsional normal modes (I, - 2,8,20). In order to determine the stability 
of the solution, six additional density models were generated, in each of which 
the period of one of the three normal modes was increased or decreased by one 
standard deviation. The superposition of the seven models is shown in Figure 4. 
It indicates that current knowledge of the torsional modes is sufficient to specify 
the density in the upper mantle but that the lower mantle and core will require 
smaller standard deviations before they can be accurately specified. 
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AN INFORMATION Til KORY APPROACH TO THE 


DENSITY OF THE EARTH 

In information theory, an important consideration in determining probability dis- 
tributions has been the maximization of missing information. 'Hiis concept is 
important because, "If we wish to avoid Introducing a bias, if we wish to avoid 
unwarranted preference for some inferences and unwarranted neglect of others, 
we must choose the compatible probability distribution with the largest value of 
the missing information," (Haierlein 1971, p. 87). The basis for this assertion 
is fully developed by Haierlein in Chapter 3 where it is shown that the missing 
information, MI, is 

Ml - £ p In p 

where p is a probability. 

David P. Rubincam (private communication 197(5) has proposed the novel idea 
that there is a relationship between the density of a body and a probability func- 
tion. 'Hie justification is developed in this manner; consider that a solid body 
is composed of infinitesimal, incompressible ping-pong balls, each having equal 
mass. Then the density of an arbitrary unit volume of the solid is proportional 
to the number of ping-pong balls in the unit volume at that point. If one of the 
ping-pong balls is painted red (all others are white), the probability of finding 
that red ping -pong ball in a given imit volume of the body is related to the number 


of ping-pong balls In that volume. Thus a relationship is developed between a 
probability function and a density function. 

Hie probability of finding the red ball is large if there are many balls in the 
volume, small if there are few, and zero if there are none. Thus the most 
general relationship between the probability and density is 

oo 

P<x.y./> = £ Uj |p(x,y,/)|' 
i * l 

For the purposes of this paper, the series is truncated after one term and it is 
asserted that 


pap 


This yields the extremal principle for the density of a solid body 


■/ 


I = / kp ln(kp) dv 


61 = 0 


( 1 ) 


with k an unspecified constant. This equation with suitable integral constraints 
forms the basis for the remainder of this paper. 


The constraints which will be applied to Equation (1) have the form of an integral 
involving the density' 


K, 




pi , (r) d\ 
pi i (r» dv 


( 2 ) 


•> 


where K, and K 2 are physical constants such as mass or moment of inertia. 
These would be incorporated Into the variational equation with Lagrange multi- 
pliers, resulting in 


/ 


J * pljtrklv + X 2 J" pi 


I ■ I kp IntkpKlv + X, I pf ( (r)dv + X 2 I pf 2 (rklv 


’Hie functional form of the density can then be determined from the Kuler- 
Lagrange equation 


dt d df 

a 0 

dp dr dp' 


where 


$ 

P 


dp 

dr 


Tliis yields 


I = 4sr 2 (kpln(kp) + X,pf,(r) + X,pf 2 (r)| 
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klnlkp) + k + Xjf|(r) + X 2 f 2 (ri = 0 



“l 

X| x 2 

p ■ 

- exp(-l ) 

* exp 1. (r) f\(r> 


k 

L k k " _ 


Hy redefining the Lagrange multipliers, 

p(r) = A expl - X, f| (r) - X 2 f 2 (r)| (3) 

’Hie number of constraints which may be incorporated is limited only by compu- 
tational considerations and is not necessarily two as in the above example. 


4 . 


3 


The values of the constants in Equation (3) are determined by substituting the 
density equation into the original constraint equations. There appears to be a 
difficulty here in that there are two equations and three unknowns, which leads 
to an Indeterminate situation. However, in practice this is overcome by always 
making the first constraint conservation of total mass 


K, 



Thus t , < r> - 1 and 

ptr) = A expl - X, - Xjf : (r)| 

Bv redefining A to inelud • exp(-X,> the equation becomes, for two constraints 

ptr) = A cxp|-X,f,(r)| 

and we are back to two equations and two unknowns. 


Rubincam's first computation of a density distribution utilized two constraints: 
mass and moment of inertia. The second constraint is placed into a spherically 
symmetric mode by noting 



♦ y-’ Klv 







+ /' ) dv 


t 


equatorial 



+ x ’ 1 * equatorial 


» 


/'■ 


^ polar + ^ oqiuiou.il 


J“ pi- ilv 0.5(( |„,| J( + - * etjiuimul * 


So that numerically 


/ 


p ilv * 5. 1 >7<» x 10* p 


gill 


j P r 1 


i- ilv * 1.5(8.051 x I0 34 ginkm*) 


with f.’(r) r- this yields the two-parameter density model 

p(r) * A expl-Xjr 2 I 

and A - 12. 14 x 1 0 1 ' gm /km ' , X, = 3. 535 x 1 0' s km'* . 


This paper is concerned with extending Rublncam's two-parameter model (un- 
published) bv developing integral constraints based on the eigenvalues u>| of 
the torsional normal modes of the earth. These constraints have the form 


w t 


j p t, (r> dv 


(I) 


Ihe precise form of f, is expressed in Appendix 4, Kquation ( A -I - 1 > . 

The function f, (r) depends on the transverse velocity v ( (r), the displacement 
potential \V(r) of the mode (see Appendix 3), and its first derivative. The re- 
quirement for v, can Ik* satisfied bv experimental seismology. Several recent 


VELOCITY (km/sec) 


papers (Choudhury, Poupinet, Perrier 1975; Braile, Keller 1975; llucstis, 
Molnar, Oliver 1975; Niazi 19731 were used to determine the mantle velocity 
profile shown In Figure 1. Qamar (1973) and Nia/I (1973) used core radii of 
3477 km and 3485 km, respectively, in 1 .lr velocity models. In Figure 1 , the 
cut-off 3481 km was chosen. This discontinuity in v, generates the core-mantle 
density discontinuities exhibited in the derived density models. 

Since f, (r) depends on the displacement potential of the normal mode, a boot- 
strap method is required. Beginning with the two-parameter density model, the 



RADIUS (1,000 km) 

Figure 1. Transverse Velocity Profile of the Mantle 


(> 
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first normal mock* Is calculated using the Hit/, method (see Appendix 3) and finite 
elements (see Appendix 2). The solution for the displacement potential is sub- 
stituted imo f, (r) and a preliminary three-pa **ameter density is generated. This 
density model generates an improved displacement potential, which in turn is 
used to generate another three-parameter density model. The procedure is con- 
tinued until the density model converges to a stable solution. No attempt was 
made to prove the uniqueness of the solution. 

Then a second normal mode was added, and the three-parameter solution was 
boot-strapped into a four-parameter solution. And this was repeated for the 
addition of a third normal mode. 


One feature of this calculation which has been temporarily ignored is the assump- 
tion used In deriving the form for the density function, Kquation (3), that f, (r) 
is independent of the density. Obviously, the displacement potential depends on 
the density model. 'Hie normalization requirement (see Appendix 3) also intro- 
tfuces a density dependence. However, within the range of density models con- 
sidered here, ihe displacement potentials varied quite slowly, indicating that the 
approximation assuming complete Independence has some value. 


The torsional normal modes chosen for this calcl ition are the zero-node oscil - 
lations of degree 2, 8, 20. The observed periods for these torsional normal 
modes are (Bullen 1975, p. 305) 


original page: IS 
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Degree 


Period (hoc) 


2 2(142. k 24. 

8 736,0 k 1.3 

20 280.8 t 0.3 

Kanamori and Anderson (1977) have shown that these observations include a 
contribution due to dispersion. Therefore, in a nun-dispersive model such as 
the present one, the values which should be used are 

Degree Period (see) 

2 2024. k 2G. 

8 727.6 t 2,8 

20 285.7* 1.1 

where the probable error has l>eon Increased to Incorporate Kanamori and 
Anderson's statement that the amplitudes of the corrections are accurate to 
20 ' . 


The numerical improvement in the boot-strap procedure can be seen from the 
calculated periods for the three torsional modes as shown in Table 1. The 
density distributions for the four models represented in Table 1 are shown in 
Figure 2. 

In order to determine the stability of the five-parameter solution, individual 
periods were Increased or decreased by the probable error and six additional 
density models were derived. These perturbed solutions with the nominal 


8 


Table 1 


Model 

Mode 

(degree) 

Calculated 

Period 

(sec) 

- ' - 

Displacement 
of Calculated 
From Observed 
Period 


2 

2535. 

-89. 

Two-Parameter 

8 



724. 3 

-3. 3 

Modes: None 

26 



286.5 

+0.8 


2 

2624. 

0 . 

Three -Parameter 

8 



747.1 

+19. 5 

Modes: 2 

26 

289.8 

+4.1 

Four-Parameter 

2 

2623. 

-1. 

Modes: 2,8 

8 

26 

727.3 

281.8 

-0.3 

-3.9 


2 

2624. 

0.0 

Five-Parameter 

8 



727. 6 

0 . 0 

Modes: 2,8,26 

26 




285.7 

0.0 


solution are shown in Figure 3 with the legend indicating which period has been 
adjusted. 

For each of the seven density models, the periods of all the zero-node torsional 
modes from L = 2 to L = 44 were calculated. These calculated periods were 
compared with the values listed by Bullen (1975) after applying a correction for 
dispersion. The calculated periods were within 1.5 standard deviations from 
the expected periods. 

The general variation of all the models is indicated in Figure 4 where all se'’en 
of the five -parameter models are superimposed. This figure indicates that 
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Figure 2. Solution Models Developed with Different 
Sets of Normal Modes 





Figure 3. Comparison of Five-Parameter Models Having Shifted Periods 

with the Nominal Solution 



Superposition of Six 
Plus the Non 



current knowledge of the periods of the torsional modes is adequate to specify 

the density distribution of the upper mantle, but not that of the lower mantle or 
core. 

Currently accepted earth density mr-aels (Hullcn, 1975) place a lower Ixmnd on 
the density at the earth's center at approximately 12 g/cc, with a probable value 
between 12 and 13 g^cc. The density then decreases as radius Increases to 
reach approximately 9.5 to 10 g^ec at the outer core boundary. A discontinuity 
is present in the core at the inner core boundary. This is not present in any of 
the models in Figure 3 because none of the constraints contains information 
about the existence of the inner core. The addition to this calculation of spheroidal 
free oscillations would include the congressional velocity, and the discontinuity 
in this velocity profile would generate an inner core discontinuity. On the mantle 
side of the core-mantle boundary, the density is believed to Ik* lK*tween 5.5 and 
(i g/cc, and then it decreases to 3.3 g/cc in the upper mantle. In Figure 3, uiree 
of the perturbed solutions are very close to this currently accepted picture of the 
earth's interior. In considering tin* differences from current earth models, it Is 
important to recall the simplicity of the present calculation and the small amount 
of experimental data included. 

lids information theory technique has taken experimentally derived numbers 
and constructed from them a nominal density model with additional perturbed 
models to indicate the general uncertainty in tiie solution. An improvement in 
the value of this technique would occur if the relative probable errors of the 
periods could be brought below 10' 4 . In order to accomplish this reduction, 
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the error :is. >* 01 * 1:1 toil with the correction for dispersion would hint* to !»«• retluccd L 

along with t lu* observational error. With the probable errors sufficiently small, 
this technique would specify very precise earth density models, and comparison 
of these with existing models would provltle a test of the assumptions u|«nn which 
this teelmique is based. 

* 1 
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APPENDIX 1 


SCALAR POTENTIALS TOR SECOND-RANK 
TANGENT TENSOR FIELDS 

A major computational simplification can bo achieved bv replacing tensor dif- 
ferential equations with equivalent scalar differential equations. Backus (19*1*5; 
19*17) suggested a technique for accomplishing this simplification in spherical 
systems. Backus' method has been followed here. More recently, Phinney and 
Burridge (Bur ridge 19*59; Phinney and Burridge 1973) and James (197*1) have 
developed similar techniques with more general applicability. 

Backus' technique is similar to replacing a vector field by three scalar potentials. 
Defining the outward normal vector n and the tangential gradient operator 

- 3 /s I d 

V = 0 — + 0 

30 sinO d$ 

then a general vector field can be written 

A = < nU(r.0,0) + V s Vtr.tf.p) - “n x V s W(r.d.0i 
in terms of three scalar potentials U, V, W. 

Backus derives a scalar representation for a general second-rank tensor In a 
three-dimensional space. However, for brevity, the result for a symmetric 
tensor i p given here, since these are of prlmarv interest. Using a tensor 


notation more in agreement with current practice. Backus' equation 4,13 (made 


symmetric) is 

T rr * P 

T lr * T r , - 0 t , ♦ e„itj k R >k 

<A1 -1) 

Tji * T tj - |l. + 8(£-H)M| K|j + 2 M.jj 
+ *inK nk \ kj ♦ fjntt " k N ki 
where i, J, k, or n = 0 or 0 

P, Q, R, L, M, N are six scalar functions of position; e |j Is the two-dimensional 
Levi-CIvita tensor; gq the metric tensor restricted to i, j = 0, <p on a unit sphere. 
R , represents tangential covariant differentiation and F (| represents two ap- 
plications of tangential covariant differentiation. 


The tangential covariant derivative is calculated in the manner of an ordinary co- 
variant derivative, except all terms introduced by the connection coefficients 
which are superscripted by "r" are eliminated and the derivatives are evaluated 
on a unit sphere so r = 1. This can be illustrated in an example by calculating 
g 1 J A l( for the basis described below, at the end of this appendix. Equation 
(Al-2) 


I d 3 

A a * A = — A ( r = I ) 

r d0 dO 

lid Id 

U = — A = — A 

r sind d0 sinfl d 0 






1-2 
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where I • 0 or 0 (not r) 


A,, ~ A l'„ A. , 



» 0 


- 0 


•V00 


1 a 2 

MI1 2 0 d0 2 


A - I 




where 1 ■ 0, or 0 (not r) 



cot0 

* - cotO 

r 



0 


A .0O 


I d- d 

A + eot0 — A 

sin*0 d0 2 dO 


K ij \ij = A , 00 + A >00 


d 2 9A I d 2 

= — A + cot# — + A 

dO 2 do sin 2 0 d0 2 


If A(r, 0,0) A(r) V v , n (0,0), where V\ m is a spherical harmonic, then 


g'i A |j = - m + I) A 


1-3 




ITils Identity hits been used In the expression for v In Equation < A 1 - 1 > where 
V(V ♦ 1)M has replaced -n 1 ' M () . For a more general equation, the reversing 
substitution may be made. 




wh0 „ c l)k are the commutation coefficient, of the boat.. 


In the bails jjlven In 


K(|u:»tlon < A1 — 2> 



for all i.j.k 


while In neneral 



* 0 




l-r> 
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APPENDIX 2 


FINITE ELEMENTS IN THE CALCULI'S OF VARIATIONS 

Hu* general format of a problem In variational calculus requires the minimization 
of I, where 

I / f(x.y,y)dx 
J '\ 

and y is a function of x over the range of integration having a first derivative y* . 
The minimization of I is carried out with respect to all functions y(x) which 
satisfy the boundary conditions specified in the problem. 

One numerical technique which can be used to put some order into the search 
for this minimum is to express y(x) as a sum of functions which satisfy the 
boundary conditions of the problem. These basis functions could be sines ard 
cosires, or if intuition gives a general picture of the solution, one could use :> 
set of functions having the general srtueture of the intuitive solution (r: mps, 
steps, etc.). Then a well-documented technique known as the Ritz method can 
be used to solve for the coefficients in this series of functions. If the basis 
functions are chosen well, the Ritz solution will be very accurate. 

The Ritz technique was followed in this paper. The basis functions used were 
finite element functions. A finite element basis function is zero over most of 
the range of the variational integral, and non-zero over a small, specific seg- 
ment (the finite element). The form of the function over this segment is arbitrary, 

2-1 
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but 1 h usually limited for numerical rensons to simple polynomials. In this work, 
linear finite elements were used. This means that the n-th basis function is de- 


fined over |a n , b n | which is a part of |x t . x > | and where 


I x 2 — x, I > I b n - a„l 


Then if c„ is in the segment |a„. b n | (not necessarily the midpoint), the n-th 
basis function is 


1',, ( x ) = 0 


x ( < x < a„ 


c n " a n 
h n ~ X 

b n " c n 


a„ < x v c„ 


c n < x < b„ 


h n < x < Xi 


At the ends of | x ( . x 2 |, the finite element is modified to 


b| - x 
b, - x, 


X| < x < h, 


h| < x < Xj 


These linear finite elements have the advantages that they are continuous, they 
are easy to handle numerically, and when summed they approximate an arbitrary 


function with a continuous function made of linear segments whose slopes change 
at each a n , c n , b n . 
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APPENDIX 3 


A VARIATIONAL FORMULATION FOR THE TORSIONAL ELASTIC 

MODES OF THE EARTH 


Backus (19fi7) has shown that for the clastic stress tensor, the six tensor poten- 
tials j?i von In Appendix 1 divide into two groups: R and N; and P, Q, L and M, 

If P, Q, L, M are non-zero, while R, N are zero, one generates spheroidal 
normal modes. If R, N are non-zero, while P, 0, L, M are zero, one gener- 
ates torsional modes. Since this paper deals with torsional modes, the general 
stress tensor will be specialized to P, Q, L, M equal to zero. 


If one postulates that the general form of the displacements of a torsional normal 
mode is 

U r = 0 
I d 


V 0 W 

si ntf >)<i> 

a 



where W is a scalar function of position, then the equations of elasticity 


c «m = T (U ».m + U m,*> 


T «m = X0 »ltm + 


m 


where 0 = gi k e 


jk 
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n-i 


(with l hr differentiation indicated by A m now general covariant differentiation! 
yield a relationship developed by Backus between tin* tensor potentials It and N 
and the displacement potential W, namely 



(lor l -'.NlusnoV I component) 



U> I) 


Phe basic equation of motion for the torsional modes of a spherically symmetric 
earth Is based iq>on the minimization of I, where 



V)d\ 


I' is the kinetic energy 



and V is the elastic energy 

V - r Vm o. 

’ , 1 Hill 

(Note that the torsional modes in a spherically symmetric earth have no gravi- 
tational contribution to the potential energy; this contribution is present for the 
spheroidal modes.) 


3-2 


Assuming that the displacement has a harmonic time dependence cxp(iu;t) then 
the integral equation of motion becomes (with W' = dW/dr) 


» /*dr |"pr 2 ^W’ 2 - WW' 

W 2 \ \ 

+ |C(t+ I) - l| — J - w 2 pr 2 W 2 J 


where V! is the degree of the normal mode. This integral formulation of the 
equation of motion can be checked by evaluating the corresponding Euler- 


Lagrange equation 


jM_ df 

aw ’ dT aw 7 


subject to the restriction that at the boundary surfaces 


The Euler-I.agrange equation yields 


, dR p 

fju)- rW + r * 3R + \2 - idk + I )| - W =0 
dr r 


which is in agreement with the sixth equation in Backus' equations 5.37 (19(57). 


'Hie boundary conditions derived with the Kuler-I.agrange equation become 


r 2 R = 0 


§§ 

§* 


Since a boundary at r = 0 for a sphere Is meaningless, this equation simplifies 
to 

R « 0 

Referring to Kquatlon < A 1 — 1 > In Appendix 1, and recalling that I* and are zero 
for torsional modes, this boundary condition Is seen to Ik* equivalent to the 
regular boundary conditions for elastic normal modes in a sphere 

l\i 3 T r0 - Tro 3 0 

Thus, solutions to the torsional normal mode problem found by minimizing l v 
will automatically satisfy the relevant boundary conditions. 

The advantage of this variational formulation is that one can use finite elements 
(sec Appendix 21 and the Hit/ method (Welnstoek 1952, Section 7-<»l to restate 
the integral variational equation as a matrix eigenvalue problem which can be 
solved extremely rapidly with modern computers. Although the differential 
equation will give the same answer, considerations of limited computer time 


dictated that the calculus of variations is the preferred mode of operation. 
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APPENDIX 4 


INTEGRAL CONSTRAINTS BASED ON TORSIONAL 

NORMAL MODES 


Additional integral constraints can be obtained for the torsional normal modes. 
In Equation (A3-1) of Appendix 3, I t contains an arbitrary constant which can 
be chosen so that the minimization procedure leaves I v ■ 0. Thus once the 
normal mode solution has been achieved one has an integral equation for the 
square of the eigenvalue 


CJ 


2 



•» 

- - WW' + |C(£ + |) - 
r 



(A 4-1 > 


with the normalization equation 



subject to the proper boundary conditions. Note that a = pvj , which gives the 
required functional form for the constraint equation (see Equation 2). Although 
this is a short-cut derivation of this equation, the end results are supported by 
VVeinstock's (1952) discussion of eigenvalues of membranes in Chapter 9 and by 
his short discussion of the similar problem in a solid in Problem 9-2. r >. The 
eigenvalue equation is put on a more general basis by the Theorem of Minimum 
Strain Energy |discussed by Sokolntkoff (195ft, p. 3H9) | since it is apparent that 
Equation (A4-1) is only 


1-1 


where V is the elastic strain energy |see Equation (A3-l)| a.ici the normalization 
equation is (where T Is the kinetic energy) 



! 


where k is an arbitrary constant chosen to make the normalization equal to one. 
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